Sampling rare fluctuations of height in the Oslo ricepile model 
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We have studied large deviations of the height of the pile from its mean value in the Oslo ricepile 
model. We sampled these very rare events with probabilities of order 1CP 100 by Monte Carlo 
simulations using importance sampling. These simulations check our qualitative arguement [Phys. 
Rev. E, 73, 021303, 2006] that in steady state of the Oslo ricepile model, the probability of large 



negative height fluctuations Ah 
held fixed, and k > 0. 



-aL about the mean varies as exp(— not L ) as L — > oo with a 



1. INTRODUCTION 

Large deviations of fluctuations in a system have been 
studied extensively 0, and recently have attracted 
much attention especially in the non-equilibrium station- 
ary states of driven systems 0, 0, 0, El- I n a recent 
paper, we have argued that in the critical slope type 
stochastic toppling models in d dimensions, the proba- 
bility of large negative height fluctuations Ah = —aL 
about the mean for a system of size L decays superexpo- 
nentially, as exp(— na 1 ^ L d+2 ) , for L — > oo, with a > 
held fixed Q. Here n > and to is an exponent defined 
by ((Ah) 2 ) ~ L 2u> . Since the arguements are plausible 
but not rigorous, it seems worthwhile to check them by 
numerical simulations. However, straight-forward sam- 
pling techniques fail in this the probabilities be- 
come very small even for fairly small L. For example, for 
L = 10, in d = 1, already the probability of the minimum 
slope configuration is of order 10~ 45 . 

In this paper, we numerically estimate the probabili- 
ties of large negative height fluctuations of the pile in the 
steady state in the one dimensional Oslo ricepile model 
using a variation of "go with the winners" strategy 
We estimate the full probability distribution function 
Pro6i(A/i) where Ah is height fluctuations of the pile 
about its mean value. This distribution has a scaling 
form L^ u) g(Ah/L u) ). The Id Oslo model is expected to 
be the Edwards- Wilkinson universality class, for which 
1/4 0. Our non-rigorous arguments in |]j imply 
that the scaling function should vary as exp(— nx A ) for 
x <C — 1 where k > 0. Our numerical data fully supports 
the theoretical expectation. 



2. DEFINITION OF THE MODEL 

The Oslo ricepile model ^J, EI is a stochastic 
sandpile-like model defined on a one dimensional lattice 
with a critical threshold value for the slope above which 
a toppling occurs, and the threshold is randomly reset 
after each toppling. Here we use an equivalent version of 
the rules as given in We consider a chain of lengh 
L. A configuration of the pile is specified by an integer 
height variable hi at each site i. The slope Zi at site i 



is defined to be hi — /ij+i, with /il+i = . Any site i 
with slope Zi = 1 is stable. Any site i with slope Zi > 3 
is said to be unstable, and relaxes by toppling. Slope 2 
can be either stable or unstable. Whenever slope at a 
site reaches the value 2 from a different value, because 
of incoming or outgoing grains, it is created as an unsta- 
ble 2 (denoted by 2). A 2 becomes stable 2 (denoted by 
2 without overbar) with probability p without any top- 
pling, or it topples with probability q = I— p. Whenever 
there is a toppling at site i, one grain is moved from the 
site i to i + If there is a toppling at the right end i = L, 
one grain goes out of the system. Grains are added only 
at site 1 and only after avalanche caused by the previous 
grain has stopped. 

The Id Oslo ricepile model has a remarkable Abelian 
property that the final height configuration does not de- 
pend on the order we topple the unstable sites Also, 
after addition of total L(L + 1) grains to any configura- 
tion, probabilities of different stable configurations are 
exactly the same as in the steady state, independent of 
the initial configuration [l2| |. The number of recurrent 
stable configurations in the critical states can be calcu- 
lated exactly, and is approximately -^( 1+ ^ ) 2L+1 for 
large L [P3j . In the steady state height profile always fluc- 
tuate between slope 1 and 2. The height h\ at the site 
1 has a stationary probability distribution, Probr,(hi), 
which is sharply peaked near its average value h\. For 
large system size L, the average height hi varies linearly 
with L, and and the fluctuations in hi scale as a sublin- 
ear power of L, with variance of hi varying as L 2u} , with 
< uo < 1. 



3. EXACT CALCULATION OF Prob^) 
SMALL L. 



FOR 



The probability distribution of hi in the steady state 
can be exactly calculated numerically for small L u sing 
the operator algebra satisfied by addition operators [12| . 
We recapitulate this briefly here. We denote any configu- 
ration by specifying slope values at all sites from i — 1 to 
i = L by a string of L integers (with or without overbar), 
e.g., |10 . . . 22). For unstable site 1 < i < L, the rules are 
as given below. 
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| . . .Zi-!,2,Z i+ l . . .} 



p\... Zi_i,2, z i+ i . . .)+q\ . . . (zj_i + 1),0, (z i+1 + 1) ■ • ■) 



I ■ ■ ■ Zj_i, 3, z i+1 ...)-> | . . . (zf_i + 1), 1, (z i+ i + 1) • ■ •) 

(!) 

with the convention that 1 = 1. At the left end, rules are 
as given above except that there is no left neighbour of 
site 1. At the right end i = L, the rules are 

. . . zz,_i, 2) -► p| . . . 2fc_i, 2) + 9| . + 1} 

|...Zi_i,3> -> |...z L _! + l,2) (2) 

Using these toppling rules repeatedly and the Abelian 
property, we can relax any unstable configuration. 

Let us now consider the state 22 . . . 22) where all the 
slopes are unstable 2's. If we add one more grain at site 
i = 1 in this state, we get the same state back (toppling 
the site with z% — 3 repeatedly) which implies that it is 
the steady state. So if we relax this configuration fully, 
we get probabilities of all the configuration in the steady 
state. For example, if we relax 1 22) for L = 2, we get the 
following sequence, 

|22) -> p|22) +<?|12) -> 




500 1000 1500 2000 2500 

L(L+l)(L+2) 

FIG. 1: (Color online) Logarithm (base 10) of the probabil- 
ity of occurrence of the minimum slope configuration, cal- 
culated exactly numerically, is plotted against system size 
L(L + 1)(L + 2) in the Oslo ricepile model. 



L 


Probability of the minimum slope 


2 


(1+P)<? 4 


3 


(l + 4p + 6p 2 + 5p 3 + 2p 4 )g 10 


4 


(1 + lOp + 45p 2 + 125p 3 + 241p 4 + 81p 9 + 18p 10 )g 20 



TABLE I: Table for the probability of the minimum slope for 
L = 2,3,4. 



p 2 \22)+pq\12)+pq\12)+q 2 \21) 
p 2 \22) + {p + p 2 )q\12) + {p + p 2 )q 2 \21) 

+ {p+p 2 )q Z \02) + (l+p)q i \ll). 

One can similarly calculate the steady state probabil- 
ities for higher L. In Table [I] we list the resulting ex- 
pression for the probability of the minimum slope config- 
uration 1 11 . . . 11) for L = 2, 3, 4. For larger L, the cal- 
culation becomes very tedious. There are two branches 
for relaxing any unstable site with z = 2 and the total 
CPU time increases as exp(L 3 ) as 0(L 3 ) relaxations of 
unstable sites are required to reach the steady state . 
We have calculated the probability of minimum slope ex- 
actly numerically for L < 12 using a simple code written 
in C. We used specific numerical values p = q = 1/2, 
to simplify the calculation, so that the probabilities are 
simple numbers and not polynomials in p. Even then, for 
L > 12, the computer time required becomes prohibitive. 

For large negative deviations, the probability becomes 
very small. Even for system size as small as L — 6 and 
L = 7, the probabilities of the minimum slope are 4.81 x 
KT 11 and 1.76 x 1CT 15 respectively. For L = 12, this 



probability is 1.23 x 1CT 55 . We argued in |7j that this 
probability is 0(q m ) with m = L(L + 1){L + 2)/6, for 
q — > 0, and is expected to decrease as exp(— n(q)L 3 ) for 
all q. In Fig.l we have plotted negative of logarithm 
of the probability of the minimum slope configuration 
versus L(L + \)(L + 2) for q = 1/2. The linear increase 
is in agreement with the theoretical expectation. 



4. MONTE CARLO SIMULATION WITH 
BIASED SAMPLING. 

In the simplest Monte Carlo algorithm to estimate 
probabilities of low-slope configurations in the steady 
state, one can start with the unstable configuration 
22... 22), and relax it by using the probabilistic top- 
pling rules until a final stable configuration is reached. 
If this is repeated many times, the fractional number of 
configurations with a given value of hi gives an estimate 
of the corresponding probability. However, this method 
is clearly unsuitable for estimating probabilities which 
are much smaller than 10 -10 . For estimating quantities 
like the probability of the minimal slope configuration in 
the steady state, this method is useless for L > 6 or so. 

Clearly, we need to implement some sort of importance 
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sampling. In the simplest implementation, one thinks of 
different possible configurations in the course of evolution 
of an avalanche at each step t as a branching tree. If we 
reach a configuration C t at a node at the i-th step, the 
probability of the process stopping is, say a(Ct). If we 
want to sample the low slope configurations, we do not 
want the process to die too soon. Then we do not select 
any of the nodes that correspond to stable configurations, 
and select one of the remaining branches with probability 
equal to their original probability, divided by the factor 
[1 — a{Ct)], and the final survival probability is estimated 
by product of such factors. 

However, this procedure also is not satisfactory for 
our problem as there are some unstable nodes all 
of whose possible resulting stable configurations have 
heights greater than the minimum height. For exam- 
ple let us consider a case for L — 5 where |21102) has 
two descendants |21102) and |21111), both stable. In the 
branching tree, these are like dead-end streets. The re- 
laxation process will eventually die if we encounter with 
such unstable configurations. But it is difficult to identify 
these directly and avoid them, without a computationally 
expensive depth-search. So the resulting process still has 
a nonzero probability of reaching such a node at the next 
step, and the overall probability of survival still decreases 
exponentially with the depth of the tree. 

We now describe an algorithm that does manage to 
avoid this problem, but at the cost of having to define 
and update one additional real variable at each site of 
the lattice. 

Algorithm for sampling rare events 

We start with a configuration with all sites unstable, 
and all Zi = 2. We imagine that there is a random num- 
ber, uniformly distributed between and 1, at each site 
i. After each toppling, the random number at the site is 
replaced by a new random number. Let the value of this 
random number at site i at the end of the update-step t 
be x(i,t). At t = 0, all x(i,0) are independent random 
variables lying between and 1. 

As it does not matter in which order we topple the 
unstable sites, we choose the following rule: at any time 
step, we first topple any unstable sites with Zi — 3, and 
reset the random number at that site. When all the un- 
stable sites are with slope 2, we topple the site having the 
largest random number, if the random number is greater 
than p. If the number is less than p, the avalanche stops. 
This constitutes the end of one update step. Note that 
one update step may involve more than one topplings at 
sites with z — 3, but there is exactly one toppling at a 
site with z — 2 in each update step. 

The problem is to determine the conditional probabil- 
ity of further evolutions for different avalanche histories, 
if we impose the condition that an avalanche does not 



stop. Note that our rule of selecting the unstable 2 for 
toppling introduces correlations between different x{i,t): 
if we know the random number at the site selected for 
toppling, random number at all other not-selected sites 
must be smaller. Our algorithm uses this correlation. We 
do not start by assigning specific values to the random 
number x{i,t) at all sites. The only information about 
x(i,t) which is known, and regularly updated, is the 
largest allowed value of x(i,t) determined by the known 
history of the avalanche. Let us call this x max (i, t). 

The toppling history can be specified by giving the site 
with z = 2 selected for toppling at each update step, and 
the random number at that site at the time of toppling. 
Let j(t) be the site selected for toppling at the t-th update 
step, and y(t) be the value of the random number at j(t) 
at the time, i.e. y(t) — x(j(t),t). We shall denote this 
sequence upto time t by % — {[j{r), y{r)] : r = 1 to t}. 
Given the sites that have been toppled, one can determine 
the set of unstable sites u(t) with z = 2, out of which 
the site with the maximum random number has to be 
selected at update-step t + 1. Since we know that x(j. t) 
for any site j £ u has not been selected for toppling 
earlier since it was reset, it must be smaller than all the 
corresponding y's selected since then. If it was reset at 
update step t prev (j, t), we must have 

(j, t) = Min{y(t') : t prev (j, t) < t' <t}. (3) 

At any update step t + 1, we first topple any sites with 
z = 3, and reset x max at the site to 1, and continue this 
till no slopes are 3. Let u(t) be the set of unstable sites 
with z = 2. It is straight forward to determine the condi- 
tional joint probability distribution of [j(t + l),y(t + 1)], 
given the condition that the avalanche does not stop, us- 
ing the information in % (see below). We then select 
one of these unstable sites in u(t) as the one having the 
largest random number, and assign the value of the ran- 
dom number at this site using the correct conditional 
probability distribution. 

The conditional probability distribution that x(j, t) is 
< x , given the value of x max (j, t), is 

x 

Prob(x(j, t) < x\x max (j, t)) = g( — — ) (4) 

%max [ Ji "J 

where g{£) as g{£) = £, for < £ < 1, and g(£) = 1 
for £ > 1. As there is no correlation between the values 
x(j,t) for different j's for the same time t, beyond that 
implied by the conditions that x(j,t) < x max (j,t), we 
must have 

Prob(y(t + l)<y\T t ) = <l> t (y) = JJ g( V -—) (5) 

If we put an additional condition that y(t + 1) > p, 
the corresponding conditional distribution of y(t + 1) is 
given by 

Prob(y(t + 1) < y\%, y(t + 1) > p) = (6) 
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In the appendix we describe the algorithm for generating 
a random number with a given distribution of the form 
Eq.©. 

Let F(t + 1) be the probability that y(t + 1) > p. 
Clearly, we have F{t + 1) equals to [1 — 3>t(p)]. i.e., 



F(t+1) = 1 



n * 

jeu(t) 



=0',*) 



(7) 



The relative weight of a particular history T t be- 
ing realized without the avalanche getting stopped is 
]~Ii'=i F(t'). We calculate the attrition factor F(t + 1) 
using Eq.Q. We then topple at the selected site j(t+ 1), 
and update the values of x max (j(t + 1)) = 1, and set 
XmaxW) = y(t + 1) for f ^ j. And repeat. 

After we start relaxing the unstable configuration 

1 222 2), the height at site 1 gradually decreases. At 

some step of relaxation, the height at first site becomes 
hi < h for the first time in the course of relaxation. We 
multiply all the previous factors, F(f)'s, upto this step 
and this product gives 



W(h) 



n f ® 



(8) 



{e t :ht>h} 



Averaging W(h) over many initial realizations, we get 
the probability of height at the first site being less than 
or equal to h, i.e., 

ProbUhi <h) = (W(h)) 

The estimate of probability of the minimum configu- 
ration is obtained by calculating the weight function 



n f ® 

{C t #|ll...ll>} 



(9) 



where the product is over all update steps t required to 
reach the minimum configuration. For different realiza- 
tions, we get different values of W m i n and, similarly as 
above, averaging over different values by taking many re- 
alizations gives us the probability of the minimum slope. 

We illustrate this procedure by a simple example. Con- 
sider a rice pile with L = 6. At t = 0, we have 
u(0) = {1,2,3,4,5,6}, as all sites are unstable. Also, 
at this stage x max = 1 for all sites. In this case, the 
probability distribution of y(l) is given by 

Prob(y(l) < y\y{l) > p) = y 6 /(l - p 6 ), for p < y < 1. 

(10) 

This can be generated as follows: select a random number 
z uniformly distributed between and 1. If z 1 / 6 > p, put 
y(l) = z 1 / 6 . If not, discard this value, and choose again. 
In this case, we get F(l) = 1 — p e . Then, we choose 
j(l) as one of the sites from u(0) at random, with equal 
probability. Say, we get = 2. Then, we topple at this 
site and assign x max = z 1 / 6 to other sites with unstable 
2's. Toppling makes sites 1 and 3 unstable, and we topple 
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FIG. 2: The frequency distribution of — log 10 (W m m) taking 
bin size = 1 and starting with 10 6 initial realizations for L = 
20. We fit the left tail of the data points with a Gaussian 



distribution with mean w 253.4 and variance a 



55. 



there as well. Whenever we topple at site with slope 3, 
we reset the x max at that site to 1. Finally, toppling at 
all unstable sites with z = 3, we get the configuration 
with u(l) = {2,3,4,5,6}, and x m ax at all these sites is 
reset to 1. So, now 

Prob(y(2) < y\y(l) > p) = y 5 /(l - p 5 ), for p < y < 1, 

(11) 

and F(2) = l- p 5 . Now we choose j(2) at random from 
u(l), say j(2) = 4. Toppling at this site induces toppling 
at other sites, and finally we get the configuration of un- 
stable sites u(2) = {1,3,4,5,6}, and we have x max = 1 
at all these sites. We now generate the variable j/(3), 
which turns out to have the same distribution as y(2). 
Then we have F(2) = 1 — p 5 . Now we choose a site from 
u(2). and so on. 



Results. 

We repeat the above procedure for many realizations 
and take the average of logarithm of the weight W(h). In 
Fig|2|we have plotted the numerically obtained distribu- 
tion of log W m in using 10 6 initial realizations for L = 20 
. It has a peak at log W m in ~ —253.4 and decays rapidly 
about the peak value. We fit the data point at the left 
of the peak value to a Gaussian distribution with mean 
fx w —253.4 and variance cr 2 w 55. It should be noted 
here that our simulation cannot accurately estimate the 
probabilities of near slope configurations for large L and 
the fractional error may be large, but the logarithm of 
the probabilities can be estimated with reasonably small 
fractional error 

As a check of our simulation algorithm, we calculated 
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FIG. 3: Exact numerical calculation and the Monte Carlo 
simulation for L < 12: Logarithm (base e) of the probabil- 
ity of the minimum slope configuration is plotted against the 
system size L in the ID Oslo ricepile model. The data is 
averaged over 10 6 realizations. 



FIG. 4: Comparision of (W m in), as computed directly from 
the simulations (blue *), and as estimated by using the Gaus- 
sian approximation in Eq |13l (light blue bullet) for different 
system sizes L. Also shown is — (log 10 W m in) (red +). 



the probability of the minimum slope configurations for 
small L and the numerical values match well with the 
values obtained from exact numerical calculation using 
the method in section. 3. For example, the probability 
of the minimum slope, for L = 5, is calculated to be 
1.475 x 10 -7 after averaging the data over 10 6 realiza- 
tions and the value is correct upto three significant dig- 
its. We have compared our results obtained from two 
procedures, i.e., the Monte Carlo simulation and exact 
numerical calculation and plotted negative of logarithm 
of the probabilities against L in Fig[3]for L < 12. 

For h near L, W(h) is a product of approximately L 3 
different factors F(t)'s, the logarithm of W(h) is a sum 
of L 3 random variables. While these variables F(t)'s 
are neither strictly independent nor they are identically 
distributed, our simulation suggest that correlations be- 
tween these factors at different times are weak, so that 
we can expect central- limit-theorem like result to hold. 
Then log [W(/i ~ L)\ may be expected to be normally 
distributed with a mean and variance both proportional 
to L 3 and W(h) would be log-normally distributed. In 
fact, the experimentally obtained probability distribution 
function of log W shows significant deviations from gaus- 
sian [ FigH]. 

Assuming that the distribution of the random variable 
X = -]n[W(h ~ L)] has the form -j±— exp[- (X 2 ~? )2 ], 
we get the probability of the minimum slope equal to 

Prob L (hi) = (W(hi)) = (er x ) (12) 

Thus, if the gaussian approximation holds for the distri- 



bution, 

2 

la(W(h)) « (\nW{h)) + y (13) 

This estimate need not be precise as the terms contribute 
significantly to (W(h ~ L)) are in the tail of the distri- 
butions and central limit theorem need not hold there. 

However our numerical estimate indicates that this ap- 
proximation is indeed very good. This is because the 
deviations from the gaussian behavior are stronger for 
smaller values of W, but these do not contribute much 
to (W). For example, from the simulation for L = 15, we 
get n w -270, a 2 w 95, ln(W mi „) « -229. The Gaus- 
sian approximation to distribution of In W m i n would give 
hx(W m in) ~ —223. In FigQ] we have compared the ac- 
tual values of logio(W m i n ) from the simulation with the 
estimate from the Gaussian approximation for different 
values of the system size L. 

In FigElwe have plotted negative of (log W m i n ) as a 
function of L and fitted it with a curve ax 3 + bx 2 + ex. 
We get a god fit for a = 0.0234 ± 0.0001, b = 0.16 ± 0.05 
and c= 0.11 ±0.2. 

Now we calculate the full probability distribution 
ProbL(hi) of height hi at site 1. We take the aver- 
age of this product over many realizations and estimate 
log[Pro6i(/ii < h)] using the Gaussian approximation. 
We thus determined ProbL{h\ — h) for different h, for 
L = 10, 20, 30, 40. The data is averaged over 10 6 initial 
realizations. In Fig|HJ we get a good scaling collapse by 
plotting 'ProbL{h\ = h) against the scaling variable 
(hi - hi)/L u where u w 1/4. Therefore Prob L {h x = h) 
a has scaling form Pro6 L (/ii) = L~ 1 / i g[(hi - (ft,i))/L 1 / 4 ] 
at the central region as well as at the tail. We see from 
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FIG. 5: Monte Carlo estimate of (logWmin) plotted against 
system size L in the ID Oslo ricepile model. The data points 



are fitted with ax 
c = 0.1. 



+ far + cx where a = 0.0234, b = 0.16 and 



FIG. 7: log 10 [Prob L (Ahi/L 1/4 < 



-x 

-4 . 



for x > 0, has been 



plotted against the scaling variable x where Ahi = (hi — hi). 
The data points are fitted with a straight line 0.11a; and it 
shows that the scaling function g{—x) varies as exp(— kx 4 ) 
for x > 1. 
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SUMMARY 

We have done a Monte Carlo simulation using impor- 
tance sampling to study large deviations in the one di- 
mensional Oslo ricepile model. We estimated probabili- 
ties of large height fluctuations of the pile and these prob- 
abilities are order 10~ 100 or even much less than this. We 
have shown that logarithm of the probability of large neg- 
ative height fluctuation Ah = —aL varies as — a 4 L 3 for 
large L. We also calculated numerically the full probabil- 
ity distribution of height of the pile, and find that it has 
scaling form L~ 1 / 4 <?(A/i/ L 1 / 4 ), with log[g(x)] varying as 
—x A for large negative x. 



FIG. 6: Scaling collapse for the probability distribution 
of height at site 1 for systems of size L = 10, 20, 30, 40. 
L 1//4 Pro6z,(A/ii) has been plotted against (hi — h^/L 1 ^ 4 . 
Each point is averaged over 10 6 realizations. 



the scaling plot the scaling function is highly asymmetric 
about the origin. 

Since the probability of minimum slope configurations 
varies as exp(— kL 3 ) for large L and the height fluctuation 
Ahi = (hi — hi) scales with system sizes as L 1 / 4 , the 
scaling function g(x) must vary as exp(— kx ). In Fig0 
we plot logarithm of Probh(Ahi/ 'L 1 / 4 < —x) versus x 4 
for x > (i.e., fluctuations below average height). The 
data gives a reasonably good fit to a straight line with 
slope 0.11. 



APPENDIX 

We illustate the procedure for generating the largest 
among x(j,tys, where j G u, using the algorithm given 
here. Consider independent random variables Xi, x<i, 
X3, X4,, X5 which are known to be uniformly distributed 
between the respective intervals as given below, 



xi e [0, 1] 

x 2 e [o,u] 

x 3 e [0,u] 

Xi G [0,v] 

x 5 g [0, v] 



where we take 1 > u > v without loss of general- 
ity. Then the cumulative probability distribution of y, 
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X ~ 



FIG. 8: The probability Prob(x max < x) of x max being less 
than x versus x. 



the largest among these random numbers, is given by 

Prob(y < x) — x, for u < x < 1 

= x 3 /u 2 , for v < x < u; 

= x 5 /(u 2 v 2 ), for < x < v. (14) 

The distribution is drawn schematically in Fig. [S] To 
generate a variable y with this distribution, we use the 
following algorithm: generate a number z randomly 
between and 1. Then following cases are possible. 

1. If u < z < 1 we choose the largest random 
number to be y = z and the maximum is surely x$. 

2. If v 3 /u 2 < z < u, we choose the largest random 
number to be y = (zu 2 ) 1 ^ 3 and the maximum is chosen 
from cc 3 , X4 and x§ with probability 1/3 each. 

3. If < z < v 3 /u 2 , we choose the largest random 
number to be y = (zv^v 2 ) 1 ^ and the maximum is 
chosen from x%, X2, £3, £4 and x$ with probability 1/5 
each. 



If the resulting value of y is less than p, the value 
is rejected, and the whole procedure is implemented 
afresh. Probability distributions for the maximum of 
more variables can be obtained similarly. 
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